1. Load packages

get packages

require(ggplot2)
require(dplyr)
require(patchwork)
require(lubridate)
require(data.table)
require(git2r)
require(stringr)
require(broom)
require(viridis)
require(car)


2. Load data

data tables

# define paths
data.dir   = "/Users/consti/Desktop/PhD/Publication_material/10_PhenogrowthEx3/Rdata/Experiment"
output.dir = "/Users/consti/Desktop/PhD/Publication_material/10_PhenogrowthEx3/Routput"

#load data
clim.data    = data.frame(fread(paste(data.dir,"climate_data.csv",sep="/"), header=T, sep=",")) # Climate 
volume.data  = data.frame(fread(paste(data.dir,"volume.csv",sep="/"), header=T, sep=","))       # Stem volume
pheno.data   = data.frame(fread(paste(data.dir,"phenology.csv",sep="/")))                       # Phenology
biomass.data = data.frame(fread(paste(data.dir,"biomass.csv",sep="/")))                         # Biomass

ggplot themes

#define plot themes
plotTheme = theme(legend.position   = "right",
                  legend.background = element_rect(fill=NA, size=0.5, linetype="solid"),
                  legend.text       = element_text(color="black"),
                  panel.grid.major  = element_blank(),
                  panel.grid.minor  = element_blank(),
                  panel.background  = element_rect(colour = NA, size=1, fill="grey97"),
                  axis.line         = element_line(color = "black"),
                  strip.background  = element_rect(fill=NA),
                  strip.text        = element_text(colour = 'black'))

plotTheme2 = theme(legend.position   = "none",
                   legend.background = element_rect(fill=NA, size=0.5, linetype="solid"),
                   legend.text       = element_text(color="black"),
                   panel.grid.major  = element_blank(),
                   panel.grid.minor  = element_blank(),
                   panel.background  = element_rect(colour = NA, size=1, fill="grey97"),
                   axis.text         = element_blank(),
                   axis.ticks        = element_blank(),
                   strip.background  = element_rect(fill=NA),
                   strip.text        = element_text(colour = 'black'))

plotTheme3 = theme(legend.position   = "none",
                  legend.background = element_rect(fill=NA, size=0.5, linetype="solid"),
                  legend.text       = element_text(color="black"),
                  panel.grid.major  = element_blank(),
                  panel.grid.minor  = element_blank(),
                  panel.background  = element_rect(colour = NA, size=1, fill="grey97"),
                  axis.line         = element_line(color = "black"),
                  strip.background  = element_rect(fill=NA),
                  strip.text        = element_text(colour = 'black'))

plotTheme4 = theme(legend.position   = "right",
                  legend.background = element_rect(fill=NA, size=0.5, linetype="solid"),
                  legend.text       = element_text(color="black"),
                  panel.grid.major  = element_blank(),
                  panel.grid.minor  = element_blank(),
                  panel.background  = element_rect(colour = NA, size=1, fill="grey97"),
                  axis.text.x       = element_text(angle = 45, hjust = 1),
                  axis.line.y       = element_line(color = "black"),
                  strip.background  = element_rect(fill=NA),
                  strip.text        = element_text(colour = 'black'))


3. Data preparation

climate data

data management

#Convert dates to DOY
clim.data$date = dmy(clim.data$date)

#get median deviations from control + interquartile ranges
clim.data1 = data.frame(clim.data %>%
                     filter(!is.na(inside_tmp)) %>% #delete summer period
                     mutate(diffOutIn  = inside_tmp-outside_tmp)) %>% #add deviation column
  group_by(treatment) %>%
  summarize(medianDiff = median(diffOutIn),
            IQrange    = IQR(diffOutIn, na.rm = FALSE, type = 7))

#get new dataframe for deviation plot
clim.data2 = data.frame(clim.data %>%
                     filter(!is.na(inside_tmp)) %>% #delete NAs
                     mutate(diffOutIn  = inside_tmp-outside_tmp)) #get deviation from control

#reshape initial dataframe
clim.data3 = melt(clim.data, id.vars = c("date","treatment"), measure.vars = c("outside_tmp","inside_tmp"))

#create treatment column
clim.data3$variable2 = ifelse(clim.data3$variable=="outside_tmp", "control", 
                       ifelse(clim.data$treatment=="autumn", "autumn", "spring"))

temperature table

data.frame(clim.data1)
##   treatment medianDiff IQrange
## 1    autumn     4.6715 1.94125
## 2    spring     4.7170 2.23475
stem volume

stem volume calculation

#order treatments
Treatment = c("spring", "autumn", "autumn-spring", "control")
volume.data$Treatment = factor(volume.data$Treatment, levels=Treatment)

#calculate volumes
volume.data$Volume = 1/3 * pi * volume.data$Height *
  ((volume.data$StemDiameter/2)^2 + 
     (volume.data$StemDiameter/2)*(volume.data$StemDiameterTop/2) + 
     volume.data$StemDiameterTop^2)
                        
#sum stem and twig volume measurements within individuals for each year
VariablesVector = c("Volume")
volume.data = data.frame(volume.data %>% 
    group_by(Species,ID,Treatment,Year) %>% 
    summarize_at(VariablesVector, sum, na.rm = TRUE))

sample sizes

table(volume.data$Treatment, volume.data$Species)/4
##                
##                 Fagus sylvatica Lonicera xylosteum Quercus robur
##   spring                      9                 10             8
##   autumn                      9                  8            10
##   autumn-spring               9                  9            10
##   control                    10                  9             9

phenology

#get RMFs
biomass.data$RSR = biomass.data$Root/biomass.data$Shoot

#delete species epithet in all tables
DeleteEpithet <- function(df) {
    df$Species = gsub(' [A-z ]*', '' , df$Species)
    return(df)
}
dfnames    = list(pheno.data,biomass.data,volume.data)
dfs        = lapply(dfnames, DeleteEpithet)
names(dfs) = c("pheno.data","biomass.data","volume.data")

#Ordering of factors
Treatments      = c("control", "autumn", "autumn-spring", "spring") #order treatments
pheno.data$Treatment = factor(pheno.data$Treatment, levels=Treatments, ordered=T)
Species         = c("Quercus", "Fagus", "Lonicera") #order species
pheno.data$Species   = factor(pheno.data$Species, levels=Species, ordered=T)

#reshape phenology table to long format
pheno.data       = melt(dfs$pheno.data, id.vars = c("Species","ID","Treatment"))

#Convert dates to DOY
pheno.data$value = dmy(pheno.data$value)
pheno.data$DOY   = yday(pheno.data$value) # Jan 1 = day 1 

#create year and phenophase columns
pheno.data$Year  = str_sub(pheno.data$variable, str_length(pheno.data$variable)-3, str_length(pheno.data$variable))
pheno.data$Pheno = str_sub(pheno.data$variable, 1, str_length(pheno.data$variable)-5)
pheno.data       = select(pheno.data, -c(variable, value))
pheno.data       = na.omit(pheno.data)

#from long (many rows) to short (multiple columns) format
pheno.data       = dcast(pheno.data, Species + ID + Treatment + Year ~ Pheno, value.var = "DOY")

#get growing season length information (leaf-out to senescence)
pheno.data$GSL   = pheno.data$Senescence - pheno.data$Leaf_out

#constrain late leafers
pheno.data$Leaf_out = ifelse(pheno.data$Leaf_out>140,140,pheno.data$Leaf_out)

biomass

#Melt data to a long format
biomass.data       = melt(dfs$biomass.data, id = c("Species", "ID", "Treatment"), 
                  variable.name="organ", value.name="biomass")

#Ordering of factors
Treatments     = c("control", "autumn", "autumn-spring", "spring") #order treatments
biomass.data$Treatment = factor(biomass.data$Treatment, levels=Treatments, ordered=T)
Species        = c("Quercus", "Fagus", "Lonicera") #order species
biomass.data$Species   = factor(biomass.data$Species, levels=Species, ordered=T)
variable       = c("Total", "Shoot", "Root", "RSR") #order organ
biomass.data$organ     = factor(biomass.data$organ, levels=variable, ordered=T)


4. Analysis

phenology

treatment effect on phenology

#################################
# Summarize annual phenology data
#################################


#get average phenology anomalies per individual
###############################################

#reshape phenology columns to long format
pheno.data = melt(pheno.data, id.vars = c("Species","ID","Treatment","Year"), 
                  variable.name="phenology", value.name="day")

#get means per species and year
phenoMean = pheno.data %>% 
    group_by(Species,phenology,Year) %>% 
    summarize(Mean = mean(day, na.rm=TRUE),
              SD = sd(day, na.rm=TRUE))

#merge data
pheno.data <-merge(pheno.data, phenoMean, all.x=T, by=c("Species","phenology","Year"))

#get change in phenology relative to mean
pheno.data = pheno.data %>% 
  mutate(phenologyAnomaly = ifelse(phenology=="Leaf_out", Mean-day, day-Mean)) %>% #add column
  dplyr::select(!c(Mean)) #delete columns

#from long (many rows) to short (multiple columns) format
pheno.data3 = dcast(pheno.data, Species + ID + Treatment + Year ~ phenology, value.var = c("phenologyAnomaly"))

#summarize years
VariablesVector = c("Leaf_out","Senescence","GSL")
pheno.data3 = data.frame(pheno.data3 %>%
    filter(!Year %in% c("2017")) %>% #delete rows based on condition
    group_by(Species,ID,Treatment) %>% 
    summarize_at(VariablesVector, mean, na.rm = TRUE))


# Add biomass information
#########################

#merge phenology and biomass tables
pheno.data3 = merge(pheno.data3, dfs$biomass.data, all.x=T, by=c("Species","ID","Treatment"))

#reshape phenology columns to long format
pheno.data3 = melt(pheno.data3, id.vars = c("Species","ID","Treatment","Total","Root","Shoot","RSR"), 
             variable.name="phenology", value.name="phenologyAnomaly")


#get average phenology dates per individual
###########################################

#from long (many rows) to short (multiple columns) format
pheno.data1 = dcast(pheno.data, Species + ID + Treatment + Year ~ phenology, value.var = c("day"))

#summarize years
VariablesVector = c("Leaf_out","Senescence","GSL")
pheno.data1 = data.frame(pheno.data1 %>%
    filter(!Year %in% c("2017")) %>% #delete rows based on condition
    group_by(Species,ID,Treatment) %>% 
    summarize_at(VariablesVector, mean, na.rm = TRUE))

#merge phenology and biomass tables
pheno.data1 = merge(pheno.data1, dfs$biomass, all.x=T, by=c("Species","ID","Treatment"))

#reshape phenology columns to long format
pheno.data1 = melt(pheno.data1, id.vars = c("Species","ID","Treatment","Total","Root","Shoot","RSR"), 
             variable.name="phenology", value.name="day")


#########################################################################################
#########################################################################################


###########################
# Multivariate linear model
###########################


#Create spring and autumn treatments
####################################

pheno.data4        = pheno.data3
pheno.data4$spring = ifelse(pheno.data4$Treatment %in% c("control","autumn"), "no","warming")
pheno.data4$autumn = ifelse(pheno.data4$Treatment %in% c("control","spring"), "no","warming")

#Multivariate linear model
##########################

resultsLM.pheno = pheno.data4 %>% 
    group_by(Species,phenology) %>% 
    do({model = lm(phenologyAnomaly ~ spring*autumn, data=.)  # create your model
    data.frame(tidy(model))}) %>%           # get coefficient info
    filter(!term %in% c("(Intercept)")) %>% # delete rows
    dplyr::select(Species, phenology, term, p.value, estimate, std.error) %>% #delete columns
    mutate_if(is.numeric, round, 2) %>%   
    mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) #add significance
#Ordering of factors
term                 = c("springwarming", "autumnwarming", "springwarming:autumnwarming") #order treatments
resultsLM.pheno$term = factor(resultsLM.pheno$term, levels=term, ordered=T)
Species                 = c("Quercus", "Fagus", "Lonicera") #order species
resultsLM.pheno$Species = factor(resultsLM.pheno$Species, levels=Species, ordered=T)
#order statistics table to match dataframe
resultsLM.pheno = resultsLM.pheno[with(resultsLM.pheno, order(resultsLM.pheno$Species, resultsLM.pheno$phenology)), ]
data.frame(resultsLM.pheno[,-c(7)])
##     Species  phenology                        term p.value estimate std.error
## 1   Quercus   Leaf_out               springwarming    0.00    16.27      1.16
## 2   Quercus   Leaf_out               autumnwarming    0.06    -2.17      1.09
## 3   Quercus   Leaf_out springwarming:autumnwarming    0.79     0.43      1.57
## 4   Quercus Senescence               springwarming    0.02   -10.77      4.56
## 5   Quercus Senescence               autumnwarming    0.00    13.67      4.31
## 6   Quercus Senescence springwarming:autumnwarming    0.39    -5.43      6.19
## 7   Quercus        GSL               springwarming    0.27     5.50      4.88
## 8   Quercus        GSL               autumnwarming    0.02    11.50      4.61
## 9   Quercus        GSL springwarming:autumnwarming    0.46    -5.00      6.63
## 10    Fagus   Leaf_out               springwarming    0.00    13.72      1.66
## 11    Fagus   Leaf_out               autumnwarming    0.00    -6.73      1.66
## 12    Fagus   Leaf_out springwarming:autumnwarming    0.71     0.89      2.37
## 13    Fagus Senescence               springwarming    0.12    -3.28      2.05
## 14    Fagus Senescence               autumnwarming    0.00    23.11      2.05
## 15    Fagus Senescence springwarming:autumnwarming    0.74    -1.00      2.94
## 16    Fagus        GSL               springwarming    0.00    10.44      2.21
## 17    Fagus        GSL               autumnwarming    0.00    16.38      2.21
## 18    Fagus        GSL springwarming:autumnwarming    0.97    -0.11      3.16
## 19 Lonicera   Leaf_out               springwarming    0.00    21.88      1.95
## 20 Lonicera   Leaf_out               autumnwarming    0.09    -3.54      2.06
## 21 Lonicera   Leaf_out springwarming:autumnwarming    0.00    -8.84      2.83
## 22 Lonicera Senescence               springwarming    0.04    -4.75      2.19
## 23 Lonicera Senescence               autumnwarming    0.00     7.69      2.32
## 24 Lonicera Senescence springwarming:autumnwarming    0.08     5.67      3.19
## 25 Lonicera        GSL               springwarming    0.00    17.13      2.77
## 26 Lonicera        GSL               autumnwarming    0.17     4.15      2.93
## 27 Lonicera        GSL springwarming:autumnwarming    0.44    -3.17      4.04
######################################################
# Check for deviation from control (one-sample T-test)
######################################################

#get means of control group
pheno.data2 = pheno.data1 %>% 
    filter(Treatment %in% c("control")) %>%
    group_by(Species,phenology) %>% 
    summarize(Mean = mean(day, na.rm=TRUE))
#merge data
pheno.data2 <-merge(pheno.data1[!pheno.data1$Treatment=="control",], 
                    pheno.data2, all.x=T, by=c("Species","phenology"))
#get change in phenology relative to control
pheno.data2$phenologyChange = pheno.data2$day - pheno.data2$Mean

#Ordering of factors
Phenology     = c("Leaf_out", "Senescence", "GSL") #order treatments
pheno.data2$phenology = factor(pheno.data2$phenology, levels=Phenology, ordered=T)
Treatments     = c("spring", "autumn", "autumn-spring") #order treatments
pheno.data2$Treatment = factor(pheno.data2$Treatment, levels=Treatments, ordered=T)
Species         = c("Quercus", "Fagus", "Lonicera") #order species
pheno.data2$Species   = factor(pheno.data2$Species, levels=Species, ordered=T)

#T-test
resultsTT = pheno.data2 %>% 
    group_by(Species,phenology,Treatment) %>% 
    do({model = t.test(.$phenologyChange) # create your model
    data.frame(tidy(model),
               tidy(shapiro.test(x = .$phenologyChange))[1,2])}) %>%           
    dplyr::select(Species, phenology, Treatment, p.value, estimate, p.value.1) %>% #delete columns
    rename(shapiroTest=p.value.1) %>%      
    mutate(significance = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) %>% #add significance
    mutate_if(is.numeric, round, 2) 
#order statistics table to match dataframe
resultsTT = resultsTT[with(resultsTT, order(resultsTT$Species, resultsTT$phenology)), ]
data.frame(resultsTT[,-c(7)])
##     Species  phenology     Treatment p.value estimate shapiroTest
## 1   Quercus   Leaf_out        spring    0.00   -16.27        0.41
## 2   Quercus   Leaf_out        autumn    0.01     2.17        0.65
## 3   Quercus   Leaf_out autumn-spring    0.00   -14.53        0.24
## 4   Quercus Senescence        spring    0.00   -10.77        0.36
## 5   Quercus Senescence        autumn    0.00    13.67        0.73
## 6   Quercus Senescence autumn-spring    0.52    -2.53        0.37
## 7   Quercus        GSL        spring    0.03     5.50        0.33
## 8   Quercus        GSL        autumn    0.01    11.50        0.69
## 9   Quercus        GSL autumn-spring    0.01    12.00        0.38
## 10    Fagus   Leaf_out        spring    0.00   -13.72        0.41
## 11    Fagus   Leaf_out        autumn    0.00     6.73        0.17
## 12    Fagus   Leaf_out autumn-spring    0.00    -7.88        0.77
## 13    Fagus Senescence        spring    0.08    -3.28        0.49
## 14    Fagus Senescence        autumn    0.00    23.11        0.02
## 15    Fagus Senescence autumn-spring    0.00    18.83        0.22
## 16    Fagus        GSL        spring    0.00    10.44        0.60
## 17    Fagus        GSL        autumn    0.00    16.38        0.31
## 18    Fagus        GSL autumn-spring    0.00    26.72        0.42
## 19 Lonicera   Leaf_out        spring    0.00   -21.88        0.08
## 20 Lonicera   Leaf_out        autumn    0.01     3.54        0.36
## 21 Lonicera   Leaf_out autumn-spring    0.00    -9.50        0.98
## 22 Lonicera Senescence        spring    0.02    -4.75        0.48
## 23 Lonicera Senescence        autumn    0.00     7.69        0.02
## 24 Lonicera Senescence autumn-spring    0.00     8.61        0.01
## 25 Lonicera        GSL        spring    0.00    17.13        0.44
## 26 Lonicera        GSL        autumn    0.00     4.15        0.36
## 27 Lonicera        GSL autumn-spring    0.00    18.11        0.98

phenology effect on total biomass

##########################
# Univariate linear models
##########################

#reshape biomass columns to long format
pheno.data3 = melt(pheno.data3, id.vars = c("Species","ID","Treatment","phenology","phenologyAnomaly"), 
             variable.name="organ", value.name="biomass")

#Transform biomass to percentage anomaly
########################################

#get means per species
pheno.data5 = pheno.data3 %>% 
    group_by(Species,organ) %>% 
    summarize(Mean = mean(biomass, na.rm=TRUE))
#merge data
pheno.data3 <-merge(pheno.data3, pheno.data5, all.x=T, by=c("Species","organ"))
#get precentage change in biomass
pheno.data3$relativeBiomass = (pheno.data3$biomass / pheno.data3$Mean -1) * 100

#order species
Species             = c("Quercus", "Fagus", "Lonicera") 
pheno.data3$Species = factor(pheno.data3$Species, levels=Species, ordered=T)

#Extract linear model info
resultsLM.pheno2 = pheno.data3 %>% 
    group_by(Species,phenology,organ) %>% 
    do({model = lm(relativeBiomass ~ phenologyAnomaly, data=.)  # create model
    data.frame(tidy(model),                # get coefficient info
               glance(model),
               tidy(shapiro.test(x = residuals(object = model)))[1,2]
               )}) %>%                     # get model info
    filter(term != "(Intercept)") %>%      # delete intercept info
    dplyr::select(Species, phenology, organ, estimate, std.error, p.value, r.squared, p.value.2) %>% # delete columns
    rename(shapiro=p.value.2) %>%          #rename columns
    mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) %>% # add column
    mutate_if(is.numeric, round, 2)

#Order phenophases
Organs            = c("Total","Shoot","Root","RSR") #order treatments
resultsLM.pheno2$organ   = factor(resultsLM.pheno2$organ, levels=Organs, ordered=T)
#order statistics table to match dataframe
resultsLM.pheno2 = resultsLM.pheno2[with(resultsLM.pheno2, 
                                         order(resultsLM.pheno2$Species, 
                                               resultsLM.pheno2$organ, 
                                               resultsLM.pheno2$phenology)), ]
data.frame(resultsLM.pheno2)
##     Species  phenology organ estimate std.error p.value r.squared shapiro sig
## 1   Quercus   Leaf_out Total     2.46      0.77    0.00      0.22    0.39  **
## 2   Quercus Senescence Total    -2.12      0.49    0.00      0.35    0.95  **
## 3   Quercus        GSL Total    -1.30      0.67    0.06      0.10    0.35   *
## 4   Quercus   Leaf_out Shoot     1.79      0.74    0.02      0.14    0.23  **
## 5   Quercus Senescence Shoot    -1.98      0.44    0.00      0.37    0.61  **
## 6   Quercus        GSL Shoot    -1.53      0.59    0.01      0.16    0.22  **
## 7   Quercus   Leaf_out  Root     3.03      0.89    0.00      0.25    0.62  **
## 8   Quercus Senescence  Root    -2.25      0.60    0.00      0.29    0.75  **
## 9   Quercus        GSL  Root    -1.10      0.81    0.18      0.05    0.02    
## 10  Quercus   Leaf_out   RSR     0.98      0.54    0.08      0.09    0.83   *
## 11  Quercus Senescence   RSR    -0.18      0.39    0.64      0.01    0.16    
## 12  Quercus        GSL   RSR     0.38      0.45    0.40      0.02    0.36    
## 13    Fagus   Leaf_out Total     0.78      0.51    0.14      0.06    0.28    
## 14    Fagus Senescence Total    -0.16      0.36    0.66      0.01    0.60    
## 15    Fagus        GSL Total     0.27      0.41    0.52      0.01    0.44    
## 16    Fagus   Leaf_out Shoot     0.42      0.52    0.43      0.02    0.31    
## 17    Fagus Senescence Shoot     0.09      0.36    0.80      0.00    0.34    
## 18    Fagus        GSL Shoot     0.37      0.40    0.36      0.02    0.50    
## 19    Fagus   Leaf_out  Root     1.24      0.60    0.05      0.11    0.18  **
## 20    Fagus Senescence  Root    -0.49      0.43    0.27      0.03    0.57    
## 21    Fagus        GSL  Root     0.13      0.50    0.80      0.00    0.07    
## 22    Fagus   Leaf_out   RSR     0.81      0.42    0.06      0.10    0.28   *
## 23    Fagus Senescence   RSR    -0.50      0.29    0.10      0.08    0.82    
## 24    Fagus        GSL   RSR    -0.14      0.34    0.68      0.00    0.94    
## 25 Lonicera   Leaf_out Total    -0.16      0.20    0.45      0.02    0.39    
## 26 Lonicera Senescence Total     0.37      0.30    0.22      0.04    0.44    
## 27 Lonicera        GSL Total     0.01      0.22    0.96      0.00    0.44    
## 28 Lonicera   Leaf_out Shoot    -0.13      0.22    0.58      0.01    0.75    
## 29 Lonicera Senescence Shoot     0.07      0.34    0.84      0.00    0.90    
## 30 Lonicera        GSL Shoot    -0.12      0.25    0.64      0.01    0.79    
## 31 Lonicera   Leaf_out  Root    -0.21      0.30    0.50      0.01    0.50    
## 32 Lonicera Senescence  Root     0.91      0.43    0.04      0.12    0.09  **
## 33 Lonicera        GSL  Root     0.24      0.33    0.48      0.02    0.18    
## 34 Lonicera   Leaf_out   RSR    -0.09      0.33    0.79      0.00    0.27    
## 35 Lonicera Senescence   RSR     0.87      0.47    0.07      0.09    0.43   *
## 36 Lonicera        GSL   RSR     0.37      0.36    0.32      0.03    0.26
#Visually inspect model assumptions 
data.assumptions = pheno.data3
data.assumptions$category = paste(data.assumptions$Species,
                                  data.assumptions$organ,
                                  data.assumptions$phenology, sep="_") #create identifier column
category.list = as.factor(unique(data.assumptions$category))       #create category vector
par(mfrow=c(2,2))                                       #set plot layout
for (category in category.list){                        #loop over categories
  tab.subset=data.assumptions[data.assumptions$category==category, ] #create table subset
  plot(lm(relativeBiomass ~ phenologyAnomaly, data=tab.subset), main=category)
}
biomass

total and per treatment biomass distributions

#Species-level plot
ggplot(data = biomass.data, mapping = aes(biomass, fill=organ)) +
  geom_histogram(bins=10, position="identity", color="black") +
  xlab("Biomass (g) [panels 1-3]                     
                                                                                           Biomass ratio [panel 4]") +
  ylab("Count") +
  scale_fill_viridis(option="E",discrete=TRUE) +
  facet_grid(Species~organ, scales="free") + 
  plotTheme4

#Treatment-level plot
ggplot(data = biomass.data[biomass.data$organ %in% c("Shoot","Root"),], mapping = aes(biomass, fill=organ)) +
  geom_histogram(bins=10, color="black") +
  xlab("Biomass (g)") +
  ylab("Count") +
  scale_fill_viridis(option="C",discrete=TRUE) +
  facet_grid(organ+Species~Treatment, scales="free") + 
  plotTheme4

treatment analysis

#######
# ANOVA
#######

#Extract model info and check ANOVA assumptions
resultsLM = biomass.data %>% 
    group_by(Species,organ) %>% 
    do({model = aov(biomass ~ Treatment, data=.)        # create your model
    data.frame(tidy(model)[1,],                         # get coefficient info
               glance(model)[1,],                       # get model info
               leveneTest=tidy(leveneTest(model))[1,4], # check Homogeneity of variances
               tidy(shapiro.test(x = residuals(object = model)))[1,2] # check for Normality
               )}) %>%           
    #select(Species, organ, p.value, r.squared, adj.r.squared, p.value.2, p.value.3) %>% #delete columns
    rename(leveneTest=p.value.2, shapiroTest=p.value.3) %>%                             #rename columns
    mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<0.1, "*", "")),
           Treatment="control")     #add significance and dummy column for plotting
data.frame(resultsLM)
##     Species organ      term df        sumsq       meansq statistic     p.value
## 1   Quercus Total Treatment  3 8.972304e+04 2.990768e+04 6.7781574 0.001097137
## 2   Quercus Shoot Treatment  3 1.352895e+04 4.509648e+03 5.2936506 0.004317898
## 3   Quercus  Root Treatment  3 3.498875e+04 1.166292e+04 6.5838558 0.001304986
## 4   Quercus   RSR Treatment  3 4.736796e-01 1.578932e-01 1.4440644 0.247684805
## 5     Fagus Total Treatment  3 1.121741e+04 3.739135e+03 1.0439249 0.386050827
## 6     Fagus Shoot Treatment  3 8.885922e+02 2.961974e+02 0.2494104 0.861165544
## 7     Fagus  Root Treatment  3 6.091089e+03 2.030363e+03 2.1809050 0.108879735
## 8     Fagus   RSR Treatment  3 1.603263e-01 5.344209e-02 1.8326934 0.160419578
## 9  Lonicera Total Treatment  3 3.100722e+02 1.033574e+02 0.2485173 0.861777929
## 10 Lonicera Shoot Treatment  3 1.906361e+02 6.354537e+01 0.3070789 0.820067194
## 11 Lonicera  Root Treatment  3 3.923083e+02 1.307694e+02 1.1726019 0.335504179
## 12 Lonicera   RSR Treatment  3 7.252300e-02 2.417433e-02 1.7831863 0.170190852
##     r.squared adj.r.squared      sigma statistic.1   p.value.1 df.1      logLik
## 1  0.38126321   0.325014412 66.4256090   6.7781574 0.001097137    4 -205.639193
## 2  0.32489040   0.263516795 29.1872856   5.2936506 0.004317898    4 -175.212262
## 3  0.37442617   0.317555822 42.0884980   6.5838558 0.001304986    4 -188.755791
## 4  0.11604443   0.035684834  0.3306652   1.4440644 0.247684805    4   -9.438122
## 5  0.08667647   0.003647058 59.8481834   1.0439249 0.386050827    4 -201.781144
## 6  0.02217097  -0.066722576 34.4614360   0.2494104 0.861165544    4 -181.358248
## 7  0.16545943   0.089592102 30.5118457   2.1809050 0.108879735    4 -176.854389
## 8  0.14281440   0.064888438  0.1707642   1.8326934 0.160419578    4   15.012317
## 9  0.02276804  -0.068847458 20.3935328   0.2485173 0.861777929    4 -157.509534
## 10 0.02798305  -0.063143535 14.3852349   0.3070789 0.820067194    4 -144.944976
## 11 0.09904344   0.014578768 10.5603384   1.1726019 0.335504179    4 -133.817484
## 12 0.14322951   0.062907271  0.1164338   1.7831863 0.170190852    4   28.453889
##          AIC       BIC     deviance df.residual leveneTest shapiroTest sig
## 1  421.27839 429.33297 1.456079e+05          33 0.40003230  0.87098652  **
## 2  360.42452 368.47911 2.811262e+04          33 0.98423338  0.99514138  **
## 3  387.51158 395.56617 5.845758e+04          33 0.04824924  0.99099527  **
## 4   28.87624  36.93083 3.608202e+00          33 0.10392883  0.66060872    
## 5  413.56229 421.61688 1.181996e+05          33 0.23193665  0.08508579    
## 6  372.71650 380.77108 3.919049e+04          33 0.14377313  0.16415740    
## 7  363.70878 371.76337 3.072210e+04          33 0.64350555  0.44789637    
## 8  -20.02463 -11.97004 9.622935e-01          33 0.49278732  0.52158902    
## 9  325.01907 332.93666 1.330868e+04          32 0.92356152  0.53347485    
## 10 299.88995 307.80755 6.621919e+03          32 0.70338321  0.69788510    
## 11 277.63497 285.55256 3.568664e+03          32 0.77713818  0.48212633    
## 12 -46.90778 -38.99018 4.338182e-01          32 0.26694805  0.66541472    
##    Treatment
## 1    control
## 2    control
## 3    control
## 4    control
## 5    control
## 6    control
## 7    control
## 8    control
## 9    control
## 10   control
## 11   control
## 12   control
# Plots to check ANOVA assumptions
biomass.data$category = paste(biomass.data$Species,biomass.data$organ, sep="_") #create identifier column
category.list = as.factor(unique(biomass.data$category))        #create category vector
par(mfrow=c(2,2))                                       #set plot layout
for (category in category.list){                        #loop over categories
  tab.subset=biomass.data[biomass.data$category==category, ]            #create table subset
  plot(aov(biomass ~ Treatment, data=tab.subset), 1, main=category) # 1. Homogeneity of variances
  plot(aov(biomass ~ Treatment, data=tab.subset), 2, main=category) # 2. Normality
}

#####
#Plot
#####

Boxp <- ggplot(data = biomass.data, mapping = aes(x = Treatment, y = biomass, fill=Treatment)) +
  geom_boxplot(position=position_dodge(0.8),coef=1e30) +
  geom_text(data = resultsLM,
            mapping = aes(x = Inf, y = Inf, hjust = 2, vjust = 2, 
                            label = paste("R2 = ", round(r.squared,2), sig, sep="")), 
              size=3.5, color="black")+
  xlab("Treatment") +
  ylab("Biomass (g)") +
  scale_fill_viridis(option="E",discrete=TRUE) +
  facet_grid(organ~Species, scales="free") + 
  plotTheme4
Boxp

Relative biomass analysis

###########################
# Multivariate linear model
###########################


#Create spring and autumn treatments
####################################

biomass.data$spring = ifelse(biomass.data$Treatment %in% c("control","autumn"), "no","warming")
biomass.data$autumn = ifelse(biomass.data$Treatment %in% c("control","spring"), "no","warming")


#Transform biomass to percentage anomaly
########################################

#get means per species
biomass.data1 = biomass.data %>% 
                group_by(Species,organ) %>% 
                summarize(Mean = mean(biomass, na.rm=TRUE))
#merge data
biomass.data1 <-merge(biomass.data, biomass.data1, all.x=T, by=c("Species","organ"))
#get precentage change in biomass
biomass.data1$relativeBiomass = (biomass.data1$biomass / biomass.data1$Mean -1) * 100


#Multivariate linear model
##########################

resultsLM.biomass = biomass.data1 %>% 
    group_by(Species,organ) %>% 
    do({model = lm(relativeBiomass ~ spring*autumn, data=.)  # create your model
    data.frame(tidy(model))}) %>%           # get coefficient info
    filter(!term %in% c("(Intercept)")) %>% # delete rows
    dplyr::select(Species, organ, term, p.value, estimate, std.error) %>% #delete columns
    mutate_if(is.numeric, round, 2) %>%   
    mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) #add significance
## `mutate_if()` ignored the following grouping variables:
## Columns `Species`, `organ`
data.frame(resultsLM.biomass)
##     Species organ                        term p.value estimate std.error sig
## 1   Quercus Total               springwarming    0.06    35.83     18.03   *
## 2   Quercus Total               autumnwarming    0.02   -40.96     17.05  **
## 3   Quercus Total springwarming:autumnwarming    0.55    14.94     24.50    
## 4   Quercus Shoot               springwarming    0.32    17.43     17.11    
## 5   Quercus Shoot               autumnwarming    0.01   -43.38     16.18  **
## 6   Quercus Shoot springwarming:autumnwarming    0.20    30.11     23.25    
## 7   Quercus  Root               springwarming    0.02    51.71     21.27  **
## 8   Quercus  Root               autumnwarming    0.06   -38.87     20.12   *
## 9   Quercus  Root springwarming:autumnwarming    0.95     1.86     28.91    
## 10  Quercus   RSR               springwarming    0.06    27.10     13.84   *
## 11  Quercus   RSR               autumnwarming    0.72     4.81     13.08    
## 12  Quercus   RSR springwarming:autumnwarming    0.23   -22.94     18.80    
## 13    Fagus Total               springwarming    0.09    20.95     12.17   *
## 14    Fagus Total               autumnwarming    0.38    10.92     12.17    
## 15    Fagus Total springwarming:autumnwarming    0.15   -25.67     17.43    
## 16    Fagus Shoot               springwarming    0.43     9.90     12.51    
## 17    Fagus Shoot               autumnwarming    0.52     8.06     12.51    
## 18    Fagus Shoot springwarming:autumnwarming    0.44   -13.94     17.92    
## 19    Fagus  Root               springwarming    0.02    35.00     14.10  **
## 20    Fagus  Root               autumnwarming    0.31    14.56     14.10    
## 21    Fagus  Root springwarming:autumnwarming    0.05   -40.59     20.19   *
## 22    Fagus   RSR               springwarming    0.03    22.66      9.90  **
## 23    Fagus   RSR               autumnwarming    0.49     6.94      9.90    
## 24    Fagus   RSR springwarming:autumnwarming    0.12   -22.38     14.19    
## 25 Lonicera Total               springwarming    0.81    -1.52      6.14    
## 26 Lonicera Total               autumnwarming    0.60     3.46      6.49    
## 27 Lonicera Total springwarming:autumnwarming    0.98     0.25      8.94    
## 28 Lonicera Shoot               springwarming    0.35    -6.47      6.77    
## 29 Lonicera Shoot               autumnwarming    0.62    -3.57      7.16    
## 30 Lonicera Shoot springwarming:autumnwarming    0.54     6.05      9.86    
## 31 Lonicera  Root               springwarming    0.42     7.26      8.83    
## 32 Lonicera  Root               autumnwarming    0.10    15.94      9.33   *
## 33 Lonicera  Root springwarming:autumnwarming    0.44   -10.06     12.85    
## 34 Lonicera   RSR               springwarming    0.14    14.36      9.37    
## 35 Lonicera   RSR               autumnwarming    0.04    21.13      9.91  **
## 36 Lonicera   RSR springwarming:autumnwarming    0.19   -18.31     13.64
########
# T-test
########

#get means of control group
###########################

biomass.data1 = biomass.data %>% 
    filter(Treatment %in% c("control")) %>%
    group_by(Species,organ) %>% 
    summarize(Mean = mean(biomass, na.rm=TRUE))
#merge data
biomass.data1 <-merge(biomass.data[!biomass.data$Treatment=="control",], biomass.data1, all.x=T, by=c("Species","organ"))
#get precentage change in biomass
biomass.data1$relativeBiomass = ifelse(biomass.data1$organ!="RSR", (biomass.data1$biomass / biomass.data1$Mean -1) * 100,
                                      (biomass.data1$biomass - biomass.data1$Mean)*100)
#Ordering of factors
Treatments     = c("spring", "autumn", "autumn-spring") #order treatments
biomass.data1$Treatment = factor(biomass.data1$Treatment, levels=Treatments, ordered=T)
Species     = c("Quercus", "Fagus", "Lonicera") #order treatments
biomass.data1$Species = factor(biomass.data1$Species, levels=Species, ordered=T)

#T-test
#######

resultsTTbiomass = biomass.data1 %>% 
    group_by(Species,organ,Treatment) %>% 
    do({model = t.test(.$relativeBiomass)        # create your model
    data.frame(tidy(model),
               tidy(shapiro.test(x = .$relativeBiomass))[1,2])}) %>%           
    select(Species, organ, Treatment, p.value, estimate, p.value.1) %>% #delete columns
    rename(shapiroTest=p.value.1) %>%   
    mutate_if(is.numeric, round, 2) %>%   
    mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) #add significance
## `mutate_if()` ignored the following grouping variables:
## Columns `Species`, `organ`, `Treatment`
#order statistics table to match dataframe
resultsTTbiomass = resultsTTbiomass[with(resultsTTbiomass, 
                                         order(resultsTTbiomass$Species, resultsTTbiomass$organ)), ]
data.frame(resultsTTbiomass)
##     Species organ     Treatment p.value estimate shapiroTest sig
## 1   Quercus Total        spring    0.08    35.60        0.59   *
## 2   Quercus Total        autumn    0.00   -40.68        0.07  **
## 3   Quercus Total autumn-spring    0.40     9.75        0.48    
## 4   Quercus Shoot        spring    0.25    16.31        0.25    
## 5   Quercus Shoot        autumn    0.01   -40.60        0.05  **
## 6   Quercus Shoot autumn-spring    0.68     3.89        0.32    
## 7   Quercus  Root        spring    0.06    54.23        0.82   *
## 8   Quercus  Root        autumn    0.00   -40.76        0.04  **
## 9   Quercus  Root autumn-spring    0.33    15.42        0.83    
## 10  Quercus   RSR        spring    0.08    31.47        0.79   *
## 11  Quercus   RSR        autumn    0.54     5.58        0.59    
## 12  Quercus   RSR autumn-spring    0.43    10.41        0.65    
## 13    Fagus Total        spring    0.04    23.08        0.17  **
## 14    Fagus Total        autumn    0.08    12.03        0.89   *
## 15    Fagus Total autumn-spring    0.65     6.83        0.07    
## 16    Fagus Shoot        spring    0.18    10.46        0.56    
## 17    Fagus Shoot        autumn    0.19     8.51        0.98    
## 18    Fagus Shoot autumn-spring    0.79     4.25        0.04    
## 19    Fagus  Root        spring    0.02    40.81        0.21  **
## 20    Fagus  Root        autumn    0.08    16.97        0.16   *
## 21    Fagus  Root autumn-spring    0.48    10.46        0.06    
## 22    Fagus   RSR        spring    0.03    17.95        0.85  **
## 23    Fagus   RSR        autumn    0.36     5.50        0.97    
## 24    Fagus   RSR autumn-spring    0.15     5.72        0.02    
## 25 Lonicera Total        spring    0.70    -1.54        0.55    
## 26 Lonicera Total        autumn    0.48     3.49        0.67    
## 27 Lonicera Total autumn-spring    0.66     2.20        0.14    
## 28 Lonicera Shoot        spring    0.14    -6.24        0.12    
## 29 Lonicera Shoot        autumn    0.56    -3.45        0.23    
## 30 Lonicera Shoot autumn-spring    0.47    -3.85        0.71    
## 31 Lonicera  Root        spring    0.31     7.96        0.52    
## 32 Lonicera  Root        autumn    0.07    17.49        0.38   *
## 33 Lonicera  Root autumn-spring    0.05    14.41        0.78   *
## 34 Lonicera   RSR        spring    0.09     8.20        0.23   *
## 35 Lonicera   RSR        autumn    0.05    12.06        0.76   *
## 36 Lonicera   RSR autumn-spring    0.01     9.81        0.72  **


5. Figures

Figure 1

#2D density plot
Dens2D = ggplot(clim.data, aes(x=outside_tmp, y=inside_tmp-outside_tmp) ) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon") +
  scale_fill_distiller(palette=4, direction=1) +
  geom_hline(aes(yintercept=4.7))+
  coord_cartesian(ylim = c(0, 10),xlim = c(-4, 19))+
   xlab("Ambient temperature (°C)") +
  ylab("Deviation from control (°C)") +
  plotTheme

#Boxplot
boxp = ggplot(clim.data2, aes(x=treatment, y=diffOutIn, fill=treatment)) + 
  geom_boxplot(outlier.shape=NA,alpha=.6) +
  coord_flip(ylim = c(0.5, 10)) +
  xlab("") +
  ylab("") +
  scale_fill_manual(values=c("orange","green4"))+
  plotTheme2

#Histogram
histo = ggplot(clim.data2, aes(x=diffOutIn, fill=treatment)) +
  geom_histogram(binwidth=.5, alpha=.6, position="identity") +
  geom_vline(data=clim.data1, aes(xintercept=medianDiff, colour=c("green4","orange3"),alpha=.8),
             linetype="dashed", size=1)+
  xlab("Deviation from control (°C)") +
  ylab("Count (number of days)") +
  scale_fill_manual(values=c("orange","green4"))+
    scale_colour_manual(values=c("orange","green3"))+
  coord_cartesian(xlim = c(0.5, 10), ylim = c(3, 60))+
  plotTheme3

#Time plot

#extract date information
phenoMean$dateMin = as.Date(phenoMean$Mean-phenoMean$SD, origin = paste0(phenoMean$Year,"-01-01"))
phenoMean$dateMax = as.Date(phenoMean$Mean+phenoMean$SD, origin = paste0(phenoMean$Year,"-01-01"))
phenoMean$Y = ifelse(phenoMean$Species=="Lonicera",30,ifelse(phenoMean$Species=="Quercus",28,26))

TempPlot = ggplot() +
  
  geom_rect(data=phenoMean[phenoMean$phenology=='Leaf_out',], 
            mapping=aes(xmin=dateMin, xmax=dateMax, ymin=-Inf, ymax=Y),
            fill=rep(c('green4','green1',"green3"),each=3), alpha=0.4) +
  
   geom_rect(data=phenoMean[phenoMean$phenology=='Senescence',], 
            mapping=aes(xmin=dateMin, xmax=dateMax, ymin=-Inf, ymax=Y),
            fill=rep(c('orange4','orange1',"orange3"),each=3), alpha=0.4) +
  
  geom_line(data=clim.data3, aes(x = date, y= value, group = variable, colour = variable2)) +
  
  geom_segment(data=phenoMean[phenoMean$phenology=='Leaf_out',], 
               aes(x = dateMin, y = Y, xend = dateMax, yend = Y),
               size=3,
               col=rep(c('green4','green1',"green3"),each=3) )+
  
  geom_segment(data=phenoMean[phenoMean$phenology=='Senescence',], 
               aes(x = dateMin, y = Y, xend = dateMax, yend = Y),
               size=3,
               col=rep(c('orange4','orange1',"orange3"),each=3) )+
  
  scale_color_manual(values=c("orange3","black","green4"))+
  coord_cartesian(ylim=c(-7,29))+
  xlab("") +
  ylab("Temperature (°C)") +
  plotTheme

#define plot layout
layout <- "
AA
AA
AA
BD
CD
CD"

#Merge plots
Fig1_ClimPlot = TempPlot + boxp + histo + Dens2D + plot_layout(design = layout)
Fig1_ClimPlot

##########
#Save PDFs
##########

pdf(paste(output.dir,"Fig1_Temperature deviation.pdf",sep="/"), width=9, height=6, useDingbats=FALSE)
Fig1_ClimPlot
dev.off()
## quartz_off_screen 
##                 2

Figure 2

#create panel labels
dat_text       = data.frame(
  label        = c("a","b","c","d","e","f","g","h","i"),
  phenology   = rep(c("Leaf_out","Senescence","GSL"),3),
  Species      = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
  term         = "springwarming")
#Plot
Fig2_LMplot = ggplot(data = resultsLM.pheno, mapping = aes(x = term, y = estimate, 
                                                fill=phenology, alpha=term)) +
  geom_bar(position=position_dodge(), stat="identity", color="black") +
    geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+
  geom_hline(yintercept=0) +
  geom_text(data = dat_text, mapping = aes(x = Inf, y = Inf, 
                                           hjust = 1, vjust = 1,
                                           label = label,
                                           fontface = "bold"))+
  stat_summary(geom = 'text', label = resultsLM.pheno$sig, 
               fun.y = mean, vjust = 2) +
  scale_alpha_discrete(range = c(1, 0.7)) +
  xlab("Phenophase") +
  ylab("Effect size (days)") +
  scale_fill_manual(values = c("green3", "orange", "blue")) +
  facet_grid(Species~phenology) + 
  plotTheme4
Fig2_LMplot

#########
#Save PDF
#########

pdf(paste(output.dir,"Fig2_PhenologyLinearTreatmentModel.pdf",sep="/"), width=7, height=7, useDingbats=FALSE)
Fig2_LMplot
dev.off()

Figure 3

volume.data$Species = factor(volume.data$Species, levels=c("Quercus robur", "Fagus sylvatica", "Lonicera xylosteum"))
pd  = position_dodge(0.2) 
Fig3_LinePlot = ggplot(volume.data, aes(x=Year, y=Volume, colour=Treatment, group=Treatment)) + 
    stat_summary(fun = mean, geom="line", position=pd) +
    stat_summary(fun.data = "mean_se", size = 0.5, position=pd) +
    labs(x = "Year", y = "Stem volume (cm3)") +
    scale_color_manual(values=c("green3", "orange", "red3", "black")) +
    coord_cartesian(ylim = c(8, 152))+
    facet_grid(~Species) +
    plotTheme 
Fig3_LinePlot

#########
#Save PDF
#########

pdf(paste(output.dir,"Fig3_LinePlot.pdf",sep="/"), width=9, height=4, useDingbats=FALSE)
Fig3_LinePlot
dev.off()

Figure 4

#Ordering of factors
term           = c("springwarming", "autumnwarming", "springwarming:autumnwarming")
resultsLM.biomass$term = factor(resultsLM.biomass$term, levels=term, ordered=T)
#create panel labels
dat_text       = data.frame(
  label        = c("a","b","c","d","e","f","g","h","i","j","k","l"),
  organ        = rep(c("Total","Shoot","Root","RSR"),3),
  Species      = c(rep("Quercus",4), rep("Fagus",4), rep("Lonicera",4)),
  term         = "springwarming")
#Plot
LMplot = ggplot(data = resultsLM.biomass, mapping = aes(x = term, y = estimate, 
                                                fill=organ, alpha=term)) +
  geom_bar(position=position_dodge(), stat="identity", color="black") +
    geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+
  geom_hline(yintercept=0) +
  geom_text(data = dat_text, mapping = aes(x = -Inf, y = -Inf, 
                                           hjust = -0.3, vjust = -1,
                                           label = label,
                                           fontface = "bold"))+
  stat_summary(geom = 'text', label = resultsLM.biomass$sig, 
               fun.y = mean, vjust = 2) +
  scale_fill_manual(values=c("orange","yellow2","red3","grey70"))+
  scale_alpha_discrete(range = c(1, 0.7)) +
  xlab("Phenophase") +
  ylab("Effect size (%)") +
  facet_grid(Species~organ) + 
  plotTheme4
LMplot

#########
#Save PDF
#########

pdf(paste(output.dir,"Fig4_BiomassLinearTreatmentModel.pdf",sep="/"), width=8, height=7, useDingbats=FALSE)
LMplot
dev.off()

Figure 5

#Effect sizes
#############

#create panel labels
dat_text      = data.frame(
  label       = c("a","b","c","d","e","f","g","h","i","j","k","l"),
  organ       = rep(c("Total","Shoot","Root","RSR"),3),
  Species     = c(rep("Quercus",4), rep("Fagus",4), rep("Lonicera",4)),
  phenology   = "Leaf_out")

#Plot
LMplot = ggplot(data = resultsLM.pheno2, mapping = aes(x = phenology, y = estimate, 
                                                fill=organ, alpha=phenology)) +
  geom_bar(position=position_dodge(), stat="identity", color="black") +
    geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+
  geom_hline(yintercept=0) +
  geom_text(data = dat_text, mapping = aes(x = -Inf, y = -Inf, 
                                           hjust = -0.1, vjust = -1,
                                           label = label,
                                           fontface = "bold"))+
  stat_summary(geom = 'text', label = resultsLM.pheno2$sig, 
               fun.y = mean, vjust = 2) +
  scale_fill_manual(values=c("orange","yellow2","red3","grey70"))+
  scale_alpha_discrete(range = c(1, 0.7)) +
  xlab("Phenophase") +
  ylab("Effect size (%/day)") +
  facet_grid(Species~organ) + 
  plotTheme4
LMplot

#Biomass
########

#order
Organs      = c("Total","Shoot","Root","RSR") #order treatments
pheno.data3$organ = factor(pheno.data3$organ, levels=Organs, ordered=T)
#create panel labels
dat_text       = data.frame(
  label        = c("a","b","c","d","e","f","g","h","i"),
  phenology    = rep(c("Leaf_out","Senescence","GSL"),3),
  Species      = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
  organ        = "Root")
BiomassPlot = ggplot(pheno.data3[pheno.data3$organ %in% c("Total","Shoot","Root"),], 
                     aes(x=phenologyAnomaly, y=relativeBiomass, colour=organ)) + 
    geom_hline(yintercept=0) +
    geom_vline(xintercept=0) +
    geom_smooth(method=lm, fullrange = T, alpha = 0.2) +
    labs(x = "Phenology anomaly (days)", y = "Biomass anomaly (%)") +
    scale_color_manual(values=c("orange","yellow2","red3"))+
    geom_text(data = dat_text, mapping = aes(x = -Inf, y = Inf, 
                                           hjust = -0.2, vjust = 1.3,
                                           label = label,
                                           fontface = "bold"), color="black")+
    geom_text(data    = resultsLM.pheno2[resultsLM.pheno2$organ=="Total",],
              mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 2.5, 
                            label = paste(round(estimate,1), " %/day", sig, sep="")), 
              size=3.5, color="orange")+
    geom_text(data    = resultsLM.pheno2[resultsLM.pheno2$organ=="Shoot",],
              mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 4.5, 
                            label = paste(round(estimate,1), " %/day", sig, sep="")), 
              size=3.5, color="yellow2")+
    geom_text(data    = resultsLM.pheno2[resultsLM.pheno2$organ=="Root",],
              mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 6.5, 
                            label = paste(round(estimate,1), " %/day", sig, sep="")), 
              size=3.5, color="red3")+
    facet_grid(Species~phenology, scales = "free") +
    plotTheme
BiomassPlot

#RSR
RSRplot = ggplot(pheno.data3[pheno.data3$organ %in% c("RSR"),], aes(x=phenologyAnomaly, y=biomass)) + 
    geom_point(size=0.2) +  
    geom_smooth(method=lm, fullrange = F) +
    labs(x = "Day", y = "Root:shoot ratio") +
    scale_color_viridis(option="viridis",discrete=TRUE) +
    geom_text(data    = resultsLM.pheno2[resultsLM.pheno2$organ=="RSR",],
              mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 2.5, 
                            label = paste("Slope = ", round(estimate,3), " day-1", sig, sep="")), 
              size=3.5, color="black")+
    facet_grid(Species~phenology, scales = "free") +
    plotTheme
RSRplot

#########
#Save PDF
#########

pdf(paste(output.dir,"Fig5_PhenologyBiomass.pdf",sep="/"), width=8, height=7, useDingbats=FALSE)
BiomassPlot
RSRplot
dev.off()

Figure S1

#create panel labels
dat_text <- data.frame(
  label       = c("a","b","c","d","e","f","g","h","i"),
  phenology   = rep(c("Leaf_out","Senescence","GSL"),3),
  Species     = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
  Treatment   = "spring")
#Plot
FigS1_PhenologyChangePlot = ggplot(data = pheno.data2, mapping = aes(x = Treatment, y = phenologyChange, 
                                                         fill=phenology, alpha=Treatment)) +
  stat_summary(fun.y = mean, 
                 geom = "bar", color="black") + 
  stat_summary(fun.data = mean_se,  
                 geom = "errorbar", width=0.2) + 
  geom_hline(yintercept=0) +
  geom_text(data = dat_text, mapping = aes(x = -Inf, y = Inf, 
                                           hjust = -0.2, vjust = 1.5,
                                           label = label,
                                           fontface = "bold"))+
  stat_summary(geom = 'text', label = resultsTT$significance, 
               fun.y = mean, vjust = 2) +
  scale_alpha_discrete(range = c(1, 0.7)) +
  xlab("Treatment") +
  ylab("Phenological change (days)") +
  scale_fill_manual(values = c("green3", "orange", "blue")) +
  facet_grid(Species~phenology) + 
  plotTheme4
FigS1_PhenologyChangePlot

#########
#Save PDF
#########

pdf(paste(output.dir,"FigS1_PhenologyChange.pdf",sep="/"), width=7, height=7, useDingbats=FALSE)
FigS1_PhenologyChangePlot
dev.off()

Figure S2

#Plot
#####

#create panel labels
dat_text <- data.frame(
  label = c("a","b","c","d","e","f","g","h","i","j","k","l"),
  organ   = rep(c("Total","Shoot","Root","RSR"),3),
  Species     = c(rep("Quercus",4), rep("Fagus",4), rep("Lonicera",4)),
  Treatment   = "spring")
#Plot
BiomassChangePlot = ggplot(data = biomass.data1, mapping = aes(x = Treatment, y = relativeBiomass, 
                                                       fill=organ, alpha=Treatment)) +
  stat_summary(fun.y = mean, 
                 geom = "bar", color="black") + 
  stat_summary(fun.data = mean_se,  
                 geom = "errorbar", width=0.2) + 
  geom_hline(yintercept=0) +
  geom_text(data = dat_text, mapping = aes(x = -Inf, y = -Inf, 
                                           hjust = -0.3, vjust = -1,
                                           label = label,
                                           fontface = "bold"))+
  stat_summary(geom = 'text', label = resultsTTbiomass$sig, 
               fun.y = max, vjust = 2.9) +
  scale_alpha_discrete(range = c(1, 0.7)) +
  coord_cartesian(ylim = c(-60, 60))+
  xlab("Treatment") +
  ylab("Biomass change (%)") +
  scale_fill_manual(values=c("orange","yellow2","red3","grey70"))+
  facet_grid(Species~organ) + 
  plotTheme4
BiomassChangePlot

#########
#Save PDF
#########

pdf(paste(output.dir,"FigS2_BiomassChange.pdf",sep="/"), width=8, height=7, useDingbats=FALSE)
BiomassChangePlot
dev.off()


6. Reproducibility

Reproducibility receipt

## datetime
Sys.time()
## [1] "2021-04-28 11:10:32 CEST"
## repository
#git2r::repository()

## sessioninfo
sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] car_3.0-6         carData_3.0-3     viridis_0.5.1     viridisLite_0.3.0
##  [5] broom_0.5.4       stringr_1.4.0     git2r_0.27.1      data.table_1.12.8
##  [9] lubridate_1.7.4   patchwork_1.0.0   dplyr_0.8.4       ggplot2_3.2.1    
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.0.0   xfun_0.12          reshape2_1.4.3     purrr_0.3.3       
##  [5] haven_2.2.0        lattice_0.20-38    colorspace_1.4-1   generics_0.0.2    
##  [9] vctrs_0.2.2        htmltools_0.4.0    yaml_2.2.1         rlang_0.4.4       
## [13] pillar_1.4.3       foreign_0.8-72     glue_1.3.1         withr_2.1.2       
## [17] RColorBrewer_1.1-2 readxl_1.3.1       plyr_1.8.5         lifecycle_0.1.0   
## [21] cellranger_1.1.0   munsell_0.5.0      gtable_0.3.0       zip_2.0.4         
## [25] evaluate_0.14      labeling_0.3       knitr_1.28         rio_0.5.16        
## [29] forcats_0.4.0      curl_4.3           Rcpp_1.0.3         scales_1.1.0      
## [33] backports_1.1.5    abind_1.4-5        farver_2.0.3       gridExtra_2.3     
## [37] hms_0.5.3          digest_0.6.23      stringi_1.4.5      openxlsx_4.1.4    
## [41] grid_3.6.2         tools_3.6.2        magrittr_1.5       lazyeval_0.2.2    
## [45] tibble_2.1.3       crayon_1.3.4       tidyr_1.0.2        pkgconfig_2.0.3   
## [49] MASS_7.3-53        assertthat_0.2.1   rmarkdown_2.1      R6_2.4.1          
## [53] nlme_3.1-142       compiler_3.6.2